# a script that calls base methylation from bismark_aa --extendedse output
# get the number of 
# at least 10 read coverage 
# at least 20 quality score per C

#use MLDBM;
#use DB_File;
use Getopt::Long;
use List::Util qw(sum);
#use Proc::ProcessTable;

#sub memory_usage {
#  my $t = new Proc::ProcessTable;
#  foreach my $got ( @{$t->table} ) {
#    next if not $got->pid eq $$;
#    return $got->size;
#  }
#}
#my $t = localtime( );warn "$t\n";warn 'memory: '. memory_usage()/1024/1024 ."\n";

# variables
my $minqual=20;
my $mincov =10;
GetOptions(
	   'minqual=i' =>\$minqual,
           'mincov=i' =>\$mincov,
	   'prefix=s' =>\$prefix,
	   'illumina=f' =>\$illumina);
sub usage
{
print STDERR <<EndOfUsage;

Usage: perl $0 [options] input_file >outFile.bed


options:
--minqual  : minquality
--mincov   : min coverage
--prefix   : prefix for output files
--illumina : illumina casava pipeline version (default:1.7), allowed options 1.5,1.6, 1.7, 1.8, 1.9
EndOfUsage
exit;
}
my $input_file=$ARGV[0];
if(@ARGV<1){usage();}

# ARRANGE OFFSET
my $offset=64;
if($illumina>1.7){$offset=33;}

# example line from bismark_aa version
#D7ZQJ5M1_0444:1:1:1318:1996#0/1	+	chr14	49602985	49603034	CNGAGAGAGAGTTCGTTTTATGTTTAGTGAGAAATGGTGTTAGGGTTTTT	CGGAGAGAGAGCCCGTCCTATGCTTAGTGAGAAATGGTGTCAGGGTCTCCAA	Z..........hxZ..hh....h.................x.....h.hh	CT	CT	]BZPPYZZUY`]^b_deeceVZcacb`^S_BBBBBBBBBBBBBBBBBBBB

# bismark output codes are:
# z unmeth CpG
# Z meth   CpG
# x unmeth CHG
# X meth   CHG
# h unmeth CHH
# H meth   CHH


#open CpG output file and print the header
# print the methylations
open (my $out,">".$prefix."_myCpG.txt");
print $out join("\t",qw(chrBase chr base strand coverage freqC freqT)),"\n";

open (my $CHHout,">".$prefix."_myCHH.txt");
print $CHHout join("\t",qw(chrBase chr base strand coverage freqC freqT nextBase)),"\n";

open (my $CHGout,">".$prefix."_myCHG.txt");
print $CHGout join("\t",qw(chrBase chr base strand coverage freqC freqT nextBase)),"\n";


#read-in file to count C bases
my %CGmethHash=(); 
my %nonCGmethHash=(); 
my %pMeth_nonCG=();
my %CHHmethHash=();
my %CHGmethHash=();

open (my $fh,$input_file);
my $lastPos  =-1;
my $lastChrom="null";
while(<$fh>)
{
    if($_=~/Bismark/){next;} # step over the header line
    chomp;
    my @cols   = split(/\t/,$_);
    my $start  = $cols[3];my $end   =$cols[4];my $strand=$cols[1];my $chr   =$cols[2];
    my @mcalls = split("",$cols[7]); # get the bismark methylation calls
    my @gbases = split("",$cols[6]);# get the genomic bases
    my @quals  = split("",$cols[10]);# get the quality scores


    #processes hashes if start-LastPos>100
    if( ($start-$lastPos>100 && $lastPos != -1 ) || ($chr ne $lastChrom && $lastChrom ne "null"  ))
    {
      processCGmethHash(\%CGmethHash,$out);
      processnonCGmethHash(\%nonCGmethHash,\%pMeth_nonCG);
      processCHmethHash(\%CHHmethHash,$CHHout);
      processCHmethHash(\%CHGmethHash,$CHGout);
      %nonCGmethHash=();
      %CGmethHash=();
      %CHHmethHash=();
      %CHGmethHash=();
    }

    # iterate over the mapped sequence
    for( my $i=0;$i< @quals; $i++) 
    {
      if ( (ord($quals[$i])-$offset) < $minqual || ( $gbases[$i] ne "C") ){next;}
      if(( $gbases[$i] eq "C" && $i==(scalar @quals)) && $gbases[$i-1].$gbases[$i].$gbases[$i+1].$gbases[$i+2] eq "CCGG" ){next;} #if last base is a C and it is a part of CCGG motif, don't call for meth
      my $key;# initilaize the hash key
      if($strand eq "+"){$key=join("|",("F",$chr,$start+$i,$gbases[$i+1]));}else{$key=join("|",("R",$chr,$end-$i,$gbases[$i+1]));}

      if($gbases[$i+1] eq "G"){ # if genomic base is CpG
	unless(exists $CGmethHash{$key}){$CGmethHash{$key}=[0,0,0];} 
	if( $mcalls[$i] eq "Z" ){$CGmethHash{$key}->[0]++;}          # update Cs
	elsif( $mcalls[$i] eq "z"){$CGmethHash{$key}->[1]++;}        # update Ts
	else{$CGmethHash{$key}->[2]++;}                              # update other bases
      }else{                    #if genomic base is non-CpG
	unless(exists $nonCGmethHash{$key}){$nonCGmethHash{$key}=[0,0,0];}

	my $isCHG;
	if(uc($mcalls[$i]) eq "X" && (! exists $CHGmethHash{$key}) ){
	  $CHGmethHash{$key}=[0,0,0];
	}
	elsif(uc($mcalls[$i]) eq "H" && (! exists $CHHmethHash{$key})){
	  $CHHmethHash{$key}=[0,0,0];
	}

	if( $mcalls[$i] eq "X" )
        {
	  $nonCGmethHash{$key}->[0]++;
	  $CHGmethHash{$key}->[0]++;
        }
	elsif($mcalls[$i] eq "H" )
	{
	  $nonCGmethHash{$key}->[0]++;
	  $CHHmethHash{$key}->[0]++;
	}
	elsif( $mcalls[$i] eq "x" ){
	  $nonCGmethHash{$key}->[1]++;
	  $CHGmethHash{$key}->[1]++;

	}
	elsif( $mcalls[$i] eq "h" ){
	   $nonCGmethHash{$key}->[1]++;
	  $CHHmethHash{$key}->[1]++;

	}
	else{
	  $nonCGmethHash{$key}->[2]++;
	  if(uc($mcalls[$i]) eq "X"){	  $CHGmethHash{$key}->[2]++;}
	  else{	  $CHHmethHash{$key}->[2]++;}
	}
      }
    }
    $lastPos=$end;
    $lastChrom=$chr;
}
close $fh;
close $out;


# get the conversion rate and write it out!!

my $numF=scalar @{$pMeth_nonCG{"F"}};
my $numR=scalar @{$pMeth_nonCG{"R"}};

if($numF==0 || $numR==0){die("\nnot enough aligments that pass 10X coverage with 20 phred square thresholds to calculate conversion rates\nEXITING....\n\n");}

my $AvFconvRate=(sum(@{$pMeth_nonCG{"F"}}))/$numF;
my $AvRconvRate=(sum(@{$pMeth_nonCG{"R"}}))/$numR;
my $AvconvRate =(sum(@{$pMeth_nonCG{"F"}})+ sum( @{$pMeth_nonCG{"R"}}))/($numF+$numR);

my @allesSchon;push @allesSchon,@{$pMeth_nonCG{"F"}},@{$pMeth_nonCG{"R"}};
my $medFconvRate=median($pMeth_nonCG{"F"});
my $medRconvRate=median($pMeth_nonCG{"R"});
my $medconvRate =median(\@allesSchon);

my $totCpG=(scalar(@allesSchon));

my $res="";
$res .= "total otherC considered (>95% C+T): $totCpG\n";
$res .= "average conversion rate = $AvconvRate\n"; 
$res .= "median conversion rate = $medconvRate\n\n"; 

$res .= "total otherC considered (Forward) (>95% C+T): $numF\n";
$res .= "average conversion rate (Forward) = $AvFconvRate\n"; 
$res .= "median conversion rate (Forward) = $medFconvRate\n\n"; 

$res .= "total otherC considered (Reverse) (>95% C+T): $numR\n";
$res .= "average conversion rate (Reverse) = $AvRconvRate\n"; 
$res .= "median conversion rate (Reverse) = $medRconvRate\n"; 

open (my $hd,">".$prefix."_conversionRate.txt");
print $hd $res;
close $hd;

###  SUBROUTINES ###################

# process a given CG methlation hash
# writes the filter passing CGs to output file
sub processCGmethHash
{
  my($CGmethHash,$out)=(@_);

  foreach my $key (keys %{$CGmethHash})
  {
        my($strand,$chr,$loc,$nextBase)=split(/\|/,$key);
	my $noCs=$CGmethHash->{$key}->[0];
	my $noTs=$CGmethHash->{$key}->[1];
	my $noOs=$CGmethHash->{$key}->[2];
	my $Cperc=sprintf("%.2f", 100*$noCs/($noTs+$noCs+$noOs) );
	my $Tperc=sprintf("%.2f", 100*$noTs/($noTs+$noCs+$noOs) );
       	if(($noTs+$noCs)/($noTs+$noCs+$noOs) > 0.9 && ($noTs+$noCs+$noOs)>=10 ){
	  print $out join("\t",($chr.".".$loc,$chr,$loc,$strand,$noCs+$noTs+$noOs,$Cperc,$Tperc)  ),"\n"; 
        }
  }
  return 1;
}

# process a given non CG methlation hash
# writes the filter passing Cs to a hash, that hash will be used to calculate conversion rate later on
sub processnonCGmethHash
{
  my($nonCGmethHash,$CTconvArray)=(@_);

  foreach my $key (keys %{$nonCGmethHash})
  {
        my($strand,$chr,$loc)=split(/\|/,$key);
	my $noCs=$nonCGmethHash->{$key}->[0];
	my $noTs=$nonCGmethHash->{$key}->[1];
	my $noOs=$nonCGmethHash->{$key}->[2];
	my $Cperc=sprintf("%.2f", 100*$noCs/($noTs+$noCs+$noOs) );
	my $Tperc=sprintf("%.2f", 100*$noTs/($noTs+$noCs+$noOs) );
       	if(($noTs+$noCs)/($noTs+$noCs+$noOs) > 0.95 && ($noTs+$noCs+$noOs)>=10 ){
	  #print join("\t",($chr.".".$loc,$chr,$loc,$strand,$noCs+$noTs+$noOs,$Cperc,$Tperc)  ),"\n"; 
	  push @{$CTconvArray->{$strand}},(($noTs*100)/($noTs+$noCs+$noOs));
        }
  }
  return 1;
}



sub processCHmethHash
{
  my($CGmethHash,$out)=(@_);

  foreach my $key (keys %{$CGmethHash})
  {
        my($strand,$chr,$loc,$nextBase)=split(/\|/,$key);
	my $noCs=$CGmethHash->{$key}->[0];
	my $noTs=$CGmethHash->{$key}->[1];
	my $noOs=$CGmethHash->{$key}->[2];
	my $Cperc=sprintf("%.2f", 100*$noCs/($noTs+$noCs+$noOs) );
	my $Tperc=sprintf("%.2f", 100*$noTs/($noTs+$noCs+$noOs) );
       	if(($noTs+$noCs)/($noTs+$noCs+$noOs) > 0.9 && ($noTs+$noCs+$noOs)>=10 ){
	  print $out join("\t",($chr.".".$loc,$chr,$loc,$strand,$noCs+$noTs+$noOs,$Cperc,$Tperc,$nextBase)  ),"\n"; 
        }
  }
  return 1;
}


# get the median value of a given array
# array of numbers
sub median {
    my $rpole = shift;
    my @pole = @$rpole;

    my $ret;

    sort(@pole);

    if( (@pole % 2) == 1 ) {
        $ret = $pole[((@pole+1) / 2)-1];
    } else {
        $ret = ($pole[(@pole / 2)-1] + $pole[@pole / 2]) / 2;
    }

    return $ret;
}

